CheckBlenAsymmetry
Xiang Ji
February 4, 2016
Load data from All.RData
rm(list=ls()) # clean up workspace
load("/Users/Xiang/GitFolders/IGCCodonSimulation/All.RData")
paml.path <- "/Users/Xiang/GitFolders/IGCCodonSimulation/"
IGC.geo.list <- c(3.0)
for (IGC.geo in IGC.geo.list){
for (localtree in 1:3){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "estimatedTau", "paml", "unrooted", "LocalTree", toString(localtree), "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", "estimatedTau", paste(toString(IGC.geo), ".0", sep = ""), "LocalTree", toString(localtree), "summary", sep = "_"), summary_mat)}
}
# Read in new PAML results
data.path <- paste(paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "estimatedTau", "paml", "unrooted", "1stTree", "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", "estimatedTau", paste(toString(IGC.geo), ".0", sep = ""), "1stTree", "summary", sep = "_"), summary_mat)
}
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "estimatedTau", "paml", "unrooted", "2ndTree", "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", "estimatedTau", paste(toString(IGC.geo), ".0", sep = ""), "2ndTree", "summary", sep = "_"), summary_mat)
}
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "10Tau", "paml", "unrooted", "1stTree", "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", "10Tau", paste(toString(IGC.geo), ".0", sep = ""), "1stTree", "summary", sep = "_"), summary_mat)
}
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "10Tau", "paml", "unrooted", "2ndTree", "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", "10Tau", paste(toString(IGC.geo), ".0", sep = ""), "2ndTree", "summary", sep = "_"), summary_mat)
}
save.image("/Users/Xiang/GitFolders/IGCCodonSimulation/All.RData")Feb 11 update
sum(PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_1_summary["ll",] < 0)## [1] 52sum(PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_2_summary["ll",] < 0)## [1] 59sum(PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_3_summary["ll",] < 0)## [1] 33hist((PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_1_summary["ll",])[PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_1_summary["ll",] < 0],
main = "lnL difference - local tree 1",
xlab = "lnL difference")
hist((PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_2_summary["ll",])[PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_2_summary["ll",] < 0],
main = "lnL difference - local tree 2",
xlab = "lnL difference")
hist((PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_3_summary["ll",])[PAML_estimatedTau_3.0_1stTree_summary["ll",] - PAML_estimatedTau_3.0_LocalTree_3_summary["ll",] < 0],
main = "lnL difference - local tree 3",
xlab = "lnL difference")
Feb 10 update
# Re-examine N4_N5, N4_mikatae
IGC.geo.list <- c(3.0, 10.0, 50.0, 100.0, 500.0)
# N4_N5
non.zero.PAML.N4.N5.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_N5", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
mean(PAML_10.0_summary["N4_N5", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
mean(PAML_50.0_summary["N4_N5", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
mean(PAML_100.0_summary["N4_N5", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
mean(PAML_500.0_summary["N4_N5", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.N5.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_N5", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
sd(PAML_10.0_summary["N4_N5", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
sd(PAML_50.0_summary["N4_N5", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
sd(PAML_100.0_summary["N4_N5", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
sd(PAML_500.0_summary["N4_N5", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.N5.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_N10", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
mean(PAML_10.0_summary["N9_N10", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
mean(PAML_50.0_summary["N9_N10", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
mean(PAML_100.0_summary["N9_N10", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
mean(PAML_500.0_summary["N9_N10", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.N5.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_N10", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
sd(PAML_10.0_summary["N9_N10", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
sd(PAML_50.0_summary["N9_N10", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
sd(PAML_100.0_summary["N9_N10", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
sd(PAML_500.0_summary["N9_N10", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
# N4_mikatae
non.zero.PAML.N4.mikatae.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_mikataeYDR418W", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
mean(PAML_10.0_summary["N4_mikataeYDR418W", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
mean(PAML_50.0_summary["N4_mikataeYDR418W", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
mean(PAML_100.0_summary["N4_mikataeYDR418W", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
mean(PAML_500.0_summary["N4_mikataeYDR418W", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.mikatae.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_mikataeYDR418W", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
sd(PAML_10.0_summary["N4_mikataeYDR418W", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
sd(PAML_50.0_summary["N4_mikataeYDR418W", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
sd(PAML_100.0_summary["N4_mikataeYDR418W", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
sd(PAML_500.0_summary["N4_mikataeYDR418W", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.mikatae.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_mikataeYEL054C", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
mean(PAML_10.0_summary["N9_mikataeYEL054C", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
mean(PAML_50.0_summary["N9_mikataeYEL054C", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
mean(PAML_100.0_summary["N9_mikataeYEL054C", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
mean(PAML_500.0_summary["N9_mikataeYEL054C", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
non.zero.PAML.N4.mikatae.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_mikataeYEL054C", (PAML_3.0_summary[4, ] > 1e-5 & PAML_3.0_summary[6,] > 1e-5)]),
sd(PAML_10.0_summary["N9_mikataeYEL054C", (PAML_10.0_summary[4, ] > 1e-5 & PAML_10.0_summary[6,] > 1e-5)]),
sd(PAML_50.0_summary["N9_mikataeYEL054C", (PAML_50.0_summary[4, ] > 1e-5 & PAML_50.0_summary[6,] > 1e-5)]),
sd(PAML_100.0_summary["N9_mikataeYEL054C", (PAML_100.0_summary[4, ] > 1e-5 & PAML_100.0_summary[6,] > 1e-5)]),
sd(PAML_500.0_summary["N9_mikataeYEL054C", (PAML_500.0_summary[4, ] > 1e-5 & PAML_500.0_summary[6,] > 1e-5)]))
# N4_N5
par(mfrow = c(2, 2))
matplot(IGC.geo.list, cbind(non.zero.PAML.N4.N5.paralog1.mean.list, non.zero.PAML.N4.N5.paralog2.mean.list, mean.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean non.zero N4.N5 estimate")
abline(h = 0.024947887926)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(non.zero.PAML.N4.N5.paralog1.sd.list, non.zero.PAML.N4.N5.paralog2.sd.list, sd.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd non.zero N4.N5 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N4_N5
matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.mean.list, PAML.N4.N5.paralog2.mean.list, mean.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.N5 estimate")
abline(h = 0.024947887926)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.sd.list, PAML.N4.N5.paralog2.sd.list, sd.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.N5 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
par(mfrow = c(2, 2))
# N4_mikatae
matplot(IGC.geo.list, cbind(non.zero.PAML.N4.mikatae.paralog1.mean.list, non.zero.PAML.N4.mikatae.paralog2.mean.list, mean.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean non.zero N4.mikatae estimate")
abline(h = 0.0566627496729)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(non.zero.PAML.N4.mikatae.paralog1.sd.list, non.zero.PAML.N4.mikatae.paralog2.sd.list, sd.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd non.zero N4.mikatae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N4_mikatae
matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.mean.list, PAML.N4.mikatae.paralog2.mean.list, mean.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.mikatae estimate")
abline(h = 0.0566627496729)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.sd.list, PAML.N4.mikatae.paralog2.sd.list, sd.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.mikatae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
Feb 5 update
Check PAML estimate of two same tree topologies but different order of taxa
# Simulation using estimated Tau value of 1.409408
# Look at difference of each estimate (Tree 1 - Tree 2)
# log likelihood
summary(PAML_estimatedTau_3.0_1stTree_summary["ll", ] - PAML_estimatedTau_3.0_2ndTree_summary["ll", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0e+00 0e+00 0e+00 1e-08 0e+00 1e-06sd(PAML_estimatedTau_3.0_1stTree_summary["ll", ] - PAML_estimatedTau_3.0_2ndTree_summary["ll", ])## [1] 9.999999e-08# kappa
summary(PAML_estimatedTau_3.0_1stTree_summary["kappa", ] - PAML_estimatedTau_3.0_2ndTree_summary["kappa", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2e-05 0e+00 0e+00 -1e-07 0e+00 1e-05sd(PAML_estimatedTau_3.0_1stTree_summary["kappa", ] - PAML_estimatedTau_3.0_2ndTree_summary["kappa", ])## [1] 5.221362e-06# omega
summary(PAML_estimatedTau_3.0_1stTree_summary["omega", ] - PAML_estimatedTau_3.0_2ndTree_summary["omega", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1e-05 0e+00 0e+00 -3e-07 0e+00 0e+00sd(PAML_estimatedTau_3.0_1stTree_summary["omega", ] - PAML_estimatedTau_3.0_2ndTree_summary["omega", ])## [1] 1.714466e-06# N0_kluyveriYDR418W
summary(PAML_estimatedTau_3.0_1stTree_summary["N0_kluyveriYDR418W", ] - PAML_estimatedTau_3.0_2ndTree_summary["N0_kluyveriYDR418W", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5e-06 0e+00 0e+00 -1e-08 0e+00 1e-06sd(PAML_estimatedTau_3.0_1stTree_summary["N0_kluyveriYDR418W", ] - PAML_estimatedTau_3.0_2ndTree_summary["N0_kluyveriYDR418W", ])## [1] 5.772628e-07# N0_N1
summary(PAML_estimatedTau_3.0_1stTree_summary["N0_N1", ] - PAML_estimatedTau_3.0_2ndTree_summary["N0_N6", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1e-06 0e+00 0e+00 -1e-08 0e+00 0e+00sd(PAML_estimatedTau_3.0_1stTree_summary["N0_N1", ] - PAML_estimatedTau_3.0_2ndTree_summary["N0_N6", ])## [1] 1e-07There seems no difference between two different tree representations
Now check PAML results on simulation with 10*Tau
# Simulation using 10 * Tau value of 1.409408 * 10 = 14.09408
# Look at difference of each estimate (Tree 1 - Tree 2)
# log likelihood
summary(PAML_10Tau_3.0_1stTree_summary["ll", ] - PAML_10Tau_3.0_2ndTree_summary["ll", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -7.7390 0.0000 0.0000 -0.2247 0.0000 5.1960sd(PAML_10Tau_3.0_1stTree_summary["ll", ] - PAML_10Tau_3.0_2ndTree_summary["ll", ])## [1] 1.69195hist(PAML_10Tau_3.0_1stTree_summary["ll", ] - PAML_10Tau_3.0_2ndTree_summary["ll", ])
# kappa
summary(PAML_10Tau_3.0_1stTree_summary["kappa", ] - PAML_10Tau_3.0_2ndTree_summary["kappa", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.99080 0.00000 0.00000 -0.06628 0.00000 0.19920sd(PAML_10Tau_3.0_1stTree_summary["kappa", ] - PAML_10Tau_3.0_2ndTree_summary["kappa", ])## [1] 0.1903289# omega
summary(PAML_10Tau_3.0_1stTree_summary["omega", ] - PAML_10Tau_3.0_2ndTree_summary["omega", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.1325000 0.0000000 0.0000000 -0.0000435 0.0000000 0.0589100sd(PAML_10Tau_3.0_1stTree_summary["omega", ] - PAML_10Tau_3.0_2ndTree_summary["omega", ])## [1] 0.01822983# N0_kluyveriYDR418W
summary(PAML_10Tau_3.0_1stTree_summary["N0_kluyveriYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N0_kluyveriYDR418W", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.18970 0.00000 0.00000 -0.02625 0.00000 0.00000sd(PAML_10Tau_3.0_1stTree_summary["N0_kluyveriYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N0_kluyveriYDR418W", ])## [1] 0.05795811# N0_N1
summary(PAML_10Tau_3.0_1stTree_summary["N0_N1", ] - PAML_10Tau_3.0_2ndTree_summary["N0_N6", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0sd(PAML_10Tau_3.0_1stTree_summary["N0_N1", ] - PAML_10Tau_3.0_2ndTree_summary["N0_N6", ])## [1] 0# N1_N2
summary(PAML_10Tau_3.0_1stTree_summary["N1_N2", ] - PAML_10Tau_3.0_2ndTree_summary["N6_N7", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.388500 0.000000 0.000000 0.003649 0.051540 0.233900sd(PAML_10Tau_3.0_1stTree_summary["N1_N2", ] - PAML_10Tau_3.0_2ndTree_summary["N6_N7", ])## [1] 0.1072574hist(PAML_10Tau_3.0_1stTree_summary["N1_N2", ] - PAML_10Tau_3.0_2ndTree_summary["N6_N7", ])
# N1_castelliiYDR418W
summary(PAML_10Tau_3.0_1stTree_summary["N1_castelliiYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N6_castelliiYDR418W", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.233900 -0.051540 0.000000 -0.006778 0.000000 0.388500sd(PAML_10Tau_3.0_1stTree_summary["N1_castelliiYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N6_castelliiYDR418W", ])## [1] 0.1068431hist(PAML_10Tau_3.0_1stTree_summary["N1_castelliiYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N6_castelliiYDR418W", ])
# N2_N3
summary(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.076140 -0.017320 0.000000 -0.007543 0.000000 0.054660sd(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])## [1] 0.02131937hist(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])
# N2_N3
summary(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.076140 -0.017320 0.000000 -0.007543 0.000000 0.054660sd(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])## [1] 0.02131937hist(PAML_10Tau_3.0_1stTree_summary["N2_N3", ] - PAML_10Tau_3.0_2ndTree_summary["N7_N8", ])
# N2_bayanusYDR418W
summary(PAML_10Tau_3.0_1stTree_summary["N2_bayanusYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N7_bayanusYDR418W", ])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.0546600 -0.0000385 0.0000000 0.0071180 0.0173400 0.0761400sd(PAML_10Tau_3.0_1stTree_summary["N2_bayanusYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N7_bayanusYDR418W", ])## [1] 0.02080385hist(PAML_10Tau_3.0_1stTree_summary["N2_bayanusYDR418W", ] - PAML_10Tau_3.0_2ndTree_summary["N7_bayanusYDR418W", ])
Now start to check branches ratio of subtree branches of paralog 1 over paralog 2 ratios
# N0_N1
target <- PAML_3.0_summary["N0_N1", ] / PAML_3.0_summary["N0_N6", ]
hist(target, main = "N0_N1"); mean(target); sd(target)
## [1] 682.5957## [1] 1311.021# N1_N2
target <- PAML_3.0_summary["N1_N2", ] / PAML_3.0_summary["N6_N7", ]
hist(target, main = "N1_N2"); mean(target); sd(target)
## [1] 1.009262## [1] 0.2219263# N1_castellii
target <- PAML_3.0_summary["N1_castelliiYDR418W", ] / PAML_3.0_summary["N6_castelliiYEL054C", ]
hist(target, main = "N1_castellii"); mean(target); sd(target)
## [1] 1.010218## [1] 0.2422185# N2_N3
target <- PAML_3.0_summary["N2_N3", ] / PAML_3.0_summary["N7_N8", ]
hist(target, main = "N2_N3"); mean(target); sd(target)
## [1] 396.4708## [1] 1982.431# N2_bayanus
target <- PAML_3.0_summary["N2_bayanusYDR418W", ] / PAML_3.0_summary["N7_bayanusYEL054C", ]
hist(target, main = "N2_bayanus"); mean(target); sd(target)
## [1] 1.135669## [1] 0.6685038# N3_N4
target <- PAML_3.0_summary["N3_N4", ] / PAML_3.0_summary["N8_N9", ]
hist(target, main = "N3_N4"); mean(target); sd(target)
## [1] 1.500644## [1] 1.810983# N3_kudriavzevii
target <- PAML_3.0_summary["N3_kudriavzeviiYDR418W", ] / PAML_3.0_summary["N8_kudriavzeviiYEL054C", ]
hist(target, main = "N3_kudriavzevii"); mean(target); sd(target)
## [1] 1.039909## [1] 0.32901# N4_N5
target <- PAML_3.0_summary["N4_N5", ] / PAML_3.0_summary["N9_N10", ]
hist(target, main = "N4_N5"); mean(target); sd(target)
## [1] 0.9775722## [1] 1.313046# N4_mikatae
target <- PAML_3.0_summary["N4_mikataeYDR418W", ] / PAML_3.0_summary["N9_mikataeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N4_mikatae"); mean(target); sd(target)
## [1] 1.866805## [1] 1.963559# N5_paradoxus
target <- PAML_3.0_summary["N5_paradoxusYDR418W", ] / PAML_3.0_summary["N10_paradoxusYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_paradoxus"); mean(target); sd(target)
## [1] 1.377084## [1] 1.532407# N5_cerevisiae
target <- PAML_3.0_summary["N5_cerevisiaeYDR418W", ] / PAML_3.0_summary["N10_cerevisiaeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_cerevisiae"); mean(target); sd(target)
## [1] 1.041587## [1] 0.5020038# N0_N1
target <- PAML_10.0_summary["N0_N1", ] / PAML_10.0_summary["N0_N6", ]
hist(target, main = "N0_N1"); mean(target); sd(target)
## [1] 808.1555## [1] 1662.984# N1_N2
target <- PAML_10.0_summary["N1_N2", ] / PAML_10.0_summary["N6_N7", ]
hist(target, main = "N1_N2"); mean(target); sd(target)
## [1] 1.071101## [1] 0.281692# N1_castellii
target <- PAML_10.0_summary["N1_castelliiYDR418W", ] / PAML_10.0_summary["N6_castelliiYEL054C", ]
hist(target, main = "N1_castellii"); mean(target); sd(target)
## [1] 1.025923## [1] 0.3345364# N2_N3
target <- PAML_10.0_summary["N2_N3", ] / PAML_10.0_summary["N7_N8", ]
hist(target, main = "N2_N3"); mean(target); sd(target)
## [1] 161.1146## [1] 1169.193# N2_bayanus
target <- PAML_10.0_summary["N2_bayanusYDR418W", ] / PAML_10.0_summary["N7_bayanusYEL054C", ]
hist(target, main = "N2_bayanus"); mean(target); sd(target)
## [1] 32.60408## [1] 313.9553# N3_N4
target <- PAML_10.0_summary["N3_N4", ] / PAML_10.0_summary["N8_N9", ]
hist(target, main = "N3_N4"); mean(target); sd(target)
## [1] 1.538383## [1] 1.689708# N3_kudriavzevii
target <- PAML_10.0_summary["N3_kudriavzeviiYDR418W", ] / PAML_10.0_summary["N8_kudriavzeviiYEL054C", ]
hist(target, main = "N3_kudriavzevii"); mean(target); sd(target)
## [1] 1.1444## [1] 0.4797683# N4_N5
target <- PAML_10.0_summary["N4_N5", ] / PAML_10.0_summary["N9_N10", ]
hist(target, main = "N4_N5"); mean(target); sd(target)
## [1] 211.0151## [1] 1481.767# N4_mikatae
target <- PAML_10.0_summary["N4_mikataeYDR418W", ] / PAML_10.0_summary["N9_mikataeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N4_mikatae"); mean(target); sd(target)
## [1] 113.2484## [1] 1108.817# N5_paradoxus
target <- PAML_10.0_summary["N5_paradoxusYDR418W", ] / PAML_10.0_summary["N10_paradoxusYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_paradoxus"); mean(target); sd(target)
## [1] 1.290973## [1] 1.112757# N5_cerevisiae
target <- PAML_10.0_summary["N5_cerevisiaeYDR418W", ] / PAML_10.0_summary["N10_cerevisiaeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_cerevisiae"); mean(target); sd(target)
## [1] 1.040702## [1] 0.568822# N0_N1
target <- PAML_50.0_summary["N0_N1", ] / PAML_50.0_summary["N0_N6", ]
hist(target, main = "N0_N1"); mean(target); sd(target)
## [1] 543.9516## [1] 1426.697# N1_N2
target <- PAML_50.0_summary["N1_N2", ] / PAML_50.0_summary["N6_N7", ]
hist(target, main = "N1_N2"); mean(target); sd(target)
## [1] 1.073097## [1] 0.2957784# N1_castellii
target <- PAML_50.0_summary["N1_castelliiYDR418W", ] / PAML_50.0_summary["N6_castelliiYEL054C", ]
hist(target, main = "N1_castellii"); mean(target); sd(target)
## [1] 1.020654## [1] 0.2728166# N2_N3
target <- PAML_50.0_summary["N2_N3", ] / PAML_50.0_summary["N7_N8", ]
hist(target, main = "N2_N3"); mean(target); sd(target)
## [1] 1.352344## [1] 1.576323# N2_bayanus
target <- PAML_50.0_summary["N2_bayanusYDR418W", ] / PAML_50.0_summary["N7_bayanusYEL054C", ]
hist(target, main = "N2_bayanus"); mean(target); sd(target)
## [1] 1.459676## [1] 1.934382# N3_N4
target <- PAML_50.0_summary["N3_N4", ] / PAML_50.0_summary["N8_N9", ]
hist(target, main = "N3_N4"); mean(target); sd(target)
## [1] 43.96184## [1] 424.0525# N3_kudriavzevii
target <- PAML_50.0_summary["N3_kudriavzeviiYDR418W", ] / PAML_50.0_summary["N8_kudriavzeviiYEL054C", ]
hist(target, main = "N3_kudriavzevii"); mean(target); sd(target)
## [1] 1.166184## [1] 0.7123563# N4_N5
target <- PAML_50.0_summary["N4_N5", ] / PAML_50.0_summary["N9_N10", ]
hist(target, main = "N4_N5"); mean(target); sd(target)
## [1] 1.407367## [1] 1.702416# N4_mikatae
target <- PAML_50.0_summary["N4_mikataeYDR418W", ] / PAML_50.0_summary["N9_mikataeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N4_mikatae"); mean(target); sd(target)
## [1] 186.1275## [1] 1843.195# N5_paradoxus
target <- PAML_50.0_summary["N5_paradoxusYDR418W", ] / PAML_50.0_summary["N10_paradoxusYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_paradoxus"); mean(target); sd(target)
## [1] 101.6888## [1] 572.1677# N5_cerevisiae
target <- PAML_50.0_summary["N5_cerevisiaeYDR418W", ] / PAML_50.0_summary["N10_cerevisiaeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_cerevisiae"); mean(target); sd(target)
## [1] 1.336023## [1] 0.9247196# N0_N1
target <- PAML_100.0_summary["N0_N1", ] / PAML_100.0_summary["N0_N6", ]
hist(target, main = "N0_N1"); mean(target); sd(target)
## [1] 702.8846## [1] 1764.081# N1_N2
target <- PAML_100.0_summary["N1_N2", ] / PAML_100.0_summary["N6_N7", ]
hist(target, main = "N1_N2"); mean(target); sd(target)
## [1] 1.049549## [1] 0.3781572# N1_castellii
target <- PAML_100.0_summary["N1_castelliiYDR418W", ] / PAML_100.0_summary["N6_castelliiYEL054C", ]
hist(target, main = "N1_castellii"); mean(target); sd(target)
## [1] 1.032278## [1] 0.2838787# N2_N3
target <- PAML_100.0_summary["N2_N3", ] / PAML_100.0_summary["N7_N8", ]
hist(target, main = "N2_N3"); mean(target); sd(target)
## [1] 89.56181## [1] 539.3573# N2_bayanus
target <- PAML_100.0_summary["N2_bayanusYDR418W", ] / PAML_100.0_summary["N7_bayanusYEL054C", ]
hist(target, main = "N2_bayanus"); mean(target); sd(target)
## [1] 1.386252## [1] 1.510943# N3_N4
target <- PAML_100.0_summary["N3_N4", ] / PAML_100.0_summary["N8_N9", ]
hist(target, main = "N3_N4"); mean(target); sd(target)
## [1] 888.6468## [1] 5070.871# N3_kudriavzevii
target <- PAML_100.0_summary["N3_kudriavzeviiYDR418W", ] / PAML_100.0_summary["N8_kudriavzeviiYEL054C", ]
hist(target, main = "N3_kudriavzevii"); mean(target); sd(target)
## [1] 1.237806## [1] 0.9738664# N4_N5
target <- PAML_100.0_summary["N4_N5", ] / PAML_100.0_summary["N9_N10", ]
hist(target, main = "N4_N5"); mean(target); sd(target)
## [1] 1.44583## [1] 3.041035# N4_mikatae
target <- PAML_100.0_summary["N4_mikataeYDR418W", ] / PAML_100.0_summary["N9_mikataeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N4_mikatae"); mean(target); sd(target)
## [1] 2.461894## [1] 3.708844# N5_paradoxus
target <- PAML_100.0_summary["N5_paradoxusYDR418W", ] / PAML_100.0_summary["N10_paradoxusYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_paradoxus"); mean(target); sd(target)
## [1] 845.0949## [1] 7179.084# N5_cerevisiae
target <- PAML_100.0_summary["N5_cerevisiaeYDR418W", ] / PAML_100.0_summary["N10_cerevisiaeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_cerevisiae"); mean(target); sd(target)
## [1] 1.472391## [1] 1.872967# N0_N1
target <- PAML_500.0_summary["N0_N1", ] / PAML_500.0_summary["N0_N6", ]
hist(target, main = "N0_N1"); mean(target); sd(target)
## [1] 390.7935## [1] 1292.949# N1_N2
target <- PAML_500.0_summary["N1_N2", ] / PAML_500.0_summary["N6_N7", ]
hist(target, main = "N1_N2"); mean(target); sd(target)
## [1] 1219.064## [1] 8571.207# N1_castellii
target <- PAML_500.0_summary["N1_castelliiYDR418W", ] / PAML_500.0_summary["N6_castelliiYEL054C", ]
hist(target, main = "N1_castellii"); mean(target); sd(target)
## [1] 1.067131## [1] 0.2984771# N2_N3
target <- PAML_500.0_summary["N2_N3", ] / PAML_500.0_summary["N7_N8", ]
hist(target, main = "N2_N3"); mean(target); sd(target)
## [1] 1309.071## [1] 6196.195# N2_bayanus
target <- PAML_500.0_summary["N2_bayanusYDR418W", ] / PAML_500.0_summary["N7_bayanusYEL054C", ]
hist(target, main = "N2_bayanus"); mean(target); sd(target)
## [1] 547.7266## [1] 4709.178# N3_N4
target <- PAML_500.0_summary["N3_N4", ] / PAML_500.0_summary["N8_N9", ]
hist(target, main = "N3_N4"); mean(target); sd(target)
## [1] 605.8042## [1] 2709.521# N3_kudriavzevii
target <- PAML_500.0_summary["N3_kudriavzeviiYDR418W", ] / PAML_500.0_summary["N8_kudriavzeviiYEL054C", ]
hist(target, main = "N3_kudriavzevii"); mean(target); sd(target)
## [1] 1.347045## [1] 1.256499# N4_N5
target <- PAML_500.0_summary["N4_N5", ] / PAML_500.0_summary["N9_N10", ]
hist(target, main = "N4_N5"); mean(target); sd(target)
## [1] 1018.097## [1] 9689.183# N4_mikatae
target <- PAML_500.0_summary["N4_mikataeYDR418W", ] / PAML_500.0_summary["N9_mikataeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N4_mikatae"); mean(target); sd(target)
## [1] 722.6354## [1] 3206.462# N5_paradoxus
target <- PAML_500.0_summary["N5_paradoxusYDR418W", ] / PAML_500.0_summary["N10_paradoxusYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_paradoxus"); mean(target); sd(target)
## [1] 140.7491## [1] 824.1837# N5_cerevisiae
target <- PAML_500.0_summary["N5_cerevisiaeYDR418W", ] / PAML_500.0_summary["N10_cerevisiaeYEL054C", ] # paralog 1 / paralog 2
hist(target, main = "N5_cerevisiae"); mean(target); sd(target)
## [1] 1.35265## [1] 1.520308